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Abstract 


In  this  paper,  we  first  describe  a  fourth  order  accurate  finite  differ¬ 
ence  discretization  for  both  the  Laplace  equation  and  the  heat  equation 
with  Dirichlet  boundary  conditions  on  irregular  domains.  In  the  case 
of  the  heat  equation  we  use  an  implicit  time  discretization  to  avoid  the 
stringent  time  step  restrictions  associated  with  explicit  schemes.  We 
then  turn  our  focus  to  the  Stefan  problem  and  construct  a  third  or¬ 
der  accurate  method  that  also  includes  an  implicit  time  discretization. 
Multidimensional  computational  results  are  presented  to  demonstrate 
the  order  accuracy  of  these  numerical  methods. 
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1  Introduction 


Various  numerical  methods  have  been  developed  to  solve  the  Stefan  prob¬ 
lem.  These  methods  need  to  be  able  to  efficiently  solve  the  heat  equation  on 
irregular  domains  and  keep  track  of  a  moving  interface  that  may  undergo 
complex  topological  changes.  The  interface  that  separates  the  two  phases 
can  be  either  explicitly  tracked  or  implicitly  captured.  The  main  disadvan¬ 
tage  of  an  explicit  approach,  e.g.  front  tracking  (see  e.g.  [16]),  is  that  special 
care  is  needed  for  topological  changes  such  as  merging  or  breaking.  More¬ 
over,  the  explicit  treatment  of  connectivity  makes  the  method  challenging 
to  extend  to  three  spatial  dimensions.  Implicit  representations  such  as  the 
level  set  method  [24,  30]  or  the  phase-field  method  [17]  embed  the  front  as  an 
isocontour  of  a  continuous  function.  Topological  changes  are  consequently 
handled  in  a  straightforward  fashion,  and  thus  these  methods  are  readily 
implemented  in  both  two  and  three  spatial  dimensions. 

Phase-field  methods  represent  the  front  implicitly  and  have  produced 
impressive  three  dimensional  results  (see  e.g.  [17,  7]).  However,  these  meth¬ 
ods  only  have  an  approximate  representation  of  the  front  location,  and  thus 
the  discretization  of  the  diffusion  field  is  less  accurate  near  the  front  resem¬ 
bling  an  enthalpy  method  [5].  Moreover,  it  is  often  challenging  to  add  new 
physics  to  the  model  since  new  asymptotic  analysis  is  often  required.  For 
more  details  on  phase  field  methods  for  the  Stefan  problem,  see  [17]  and  the 
references  therein. 

In  this  paper,  we  employ  the  sharp  interface  implicit  representation  of 
the  level  set  method  [24,  30].  The  earliest  level  set  method  for  solidification 
type  problems  was  presented  in  [31]  where  the  authors  recast  the  equations 
of  motion  into  a  boundary  integral  equation  and  used  the  level  set  method 
to  update  the  location  of  the  interface.  In  [2] ,  the  boundary  integral  equa¬ 
tions  were  avoided  by  using  a  finite  difference  method  to  solve  the  diffusion 
equation  on  a  Cartesian  grid  with  Dirichlet  boundary  conditions  imposed 
on  the  interface.  The  jump  in  the  first  derivatives  of  the  temperature  was 
used  to  compute  an  interface  velocity  that  was  extended  to  a  band  about  the 
interface  and  used  to  evolve  the  level  set  function  in  time.  [10]  showed  that 
this  discretization  produces  results  in  accordance  with  solvability  theory.  In 
[35],  the  authors  discretized  the  heat  equation  on  a  Cartesian  grid  in  a  man¬ 
ner  quite  similar  to  that  proposed  in  [2]  resulting  in  a  nonsymmetric  linear 
system  when  applying  implicit  time  discretization.  [35]  used  front  tracking 
to  update  the  location  of  the  interface  improving  upon  the  front  tracking 
approach  proposed  in  [16]  which  used  the  smeared  out  immersed  boundary 
method  [25]  and  an  explicit  time  discretization. 
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In  [15],  the  authors  solved  a  variable  coefficient  Poisson  equation  in  the 
presence  of  an  irregular  interface  where  Dirichlet  boundary  conditions  were 
imposed.  They  used  a  finite  volume  method  that  results  in  a  nonsymmetric 
discretization  matrix.  Both  multigrid  methods  and  adaptive  mesh  refine¬ 
ment  were  used,  and  in  [14]  this  nonsymmetric  discretization  was  coupled 
to  a  volume  of  fluid  front  tracking  method  in  order  to  solve  the  Stefan  prob¬ 
lem.  In  [27]  the  authors  used  adaptive  finite  element  methods  for  both  the 
heat  equation  and  for  the  interface  evolution  producing  stunning  three  di¬ 
mensional  results.  Other  remarkable  three  dimensional  results  can  be  found 
in  the  finite  difference  diffusion  Monte  Carlo  method  of  [26].  Recently,  [36] 
formulated  a  second  order  accurate  method  for  the  Stefan  problem  in  two 
spatial  dimensions  using  a  Galerkin  finite  element  approach  to  solve  for  the 
energy  equation.  In  this  work,  the  interface  was  tracked  with  a  set  of  marker 
particles  making  the  method  potentially  hard  to  extend  to  three  spatial  di¬ 
mensions.  Moreover,  the  velocity  was  computed  under  the  assumption  that 
the  interface  cuts  the  element  in  a  straight  line.  The  interested  reader  is 
referred  to  [16],  [2],  [8]  and  the  references  therein  for  an  extensive  summary 
of  computational  results  for  the  Stefan  problem. 

Standard  convergence  proofs  use  stability  and  consistency  analysis  to 
imply  convergence,  i.e.  given  stability,  a  sufficient  condition  for  a  scheme  to 
be  pth  order  accurate  is  that  the  local  truncation  error  is  pth  order.  However, 
a  pth  order  local  truncation  error  is  not  a  necessary  condition  and  one  can 
derive  pth  order  accurate  schemes  despite  the  fact  that  their  local  truncation 
error  is  of  lower  order.  [23]  (see  also  [19])  made  this  point  in  the  context 
of  second  order,  scalar  boundary  value  problems  on  nonuniform  meshes.  In 
fact,  in  the  process  of  constructing  second  order  accurate  methods  for  such 
problems,  many  authors  had  unnecessarily  focused  on  imposing  special  re¬ 
strictions  on  the  mesh  size  in  order  to  obtain  a  second  order  accurate  local 
truncation  error,  see  e.g  [4,  11].  In  the  case  of  our  present  work,  authors  have 
also  been  misled  by  the  limitation  of  standard  convergence  analysis  proofs 
and  have  proposed  unnecessarily  complex  schemes.  For  example,  in  [2]  (see 
also  [10])  the  authors  approximate  the  Laplace  operator  with  the  standard 
second  order  accurate  central  scheme  limiting  the  overall  solution  to  second 
order  accuracy.  However,  their  discretization  of  the  Laplace  operator  for 
grid  nodes  neighboring  the  interface  amounts  to  differentiating  a  quadratic 
interpolant  of  the  temperature  twice  in  each  spatial  dimension.  [9]  refor¬ 
mulated  the  interface  treatment  with  the  use  of  ghost  cells  (based  on  the 
ghost  fluid  method  [6])  defined  by  extrapolation  of  the  temperature  across 
the  interface,  and  showed  that  local  linear  extrapolation  is  enough  to  obtain 
second  order  spatial  accuracy  for  both  the  Laplace  and  heat  equations  on 
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irregular  domains.  Moreover,  their  discretization  had  the  added  benefit  of 
yielding  a  symmetric  linear  system  as  opposed  to  the  nonsymmetric  system 
in  [2],  This  scheme  served  as  the  basis  of  a  simple  method  to  solve  the  Stefan 
problem.  It  was  further  used  in  [8]  to  show  that  one  could  obtain  solutions  in 
agreement  with  solvability  theory,  and  could  simulate  many  of  the  physical 
features  of  crystal  growth  such  as  molecular  kinetics  and  surface  tension. 

In  this  paper,  we  exploit  the  methodology  of  [9]  to  derive  a  fourth  order 
accurate  finite  difference  discretization  for  the  Laplace  equation  on  irregular 
domains.  Then  we  apply  this  framework  to  derive  a  fourth  order  accurate 
discretization  for  the  heat  equation  with  Dirichlet  boundary  conditions  on 
arbitrary  domains.  In  this  case  we  use  an  implicit  time  discretization  to 
avoid  the  stringent  time  step  restrictions  induced  by  explicit  schemes.  We 
then  turn  our  focus  to  the  Stefan  problem  with  Dirichlet  boundary  con¬ 
ditions  and  construct  a  third  order  accurate  discretization  that  includes 
implicit  integration  in  time.  Multidimensional  computational  results  are 
presented  to  verify  the  order  accuracy  of  these  numerical  methods. 

2  Laplace  Equation 

Consider  a  Cartesian  computational  domain,  P  C  R" ,  with  exterior  bound¬ 
ary,  90,  and  a  lower  dimensional  interface,  T,  that  divides  the  computational 
domain  into  disjoint  pieces,  O”  and  0+.  The  Laplace  equation  is  given  by 

A T(x)  =  f  (x) ,  x  €  O",  (1) 

where  x  =  (x,  y,  z )  is  the  vector  of  spatial  coordinates,  A  =  ^  is 

the  Laplace  operator,  and  T  is  assumed  to  be  smooth  on  P-.  On  T,  Dirichlet 
boundary  conditions  are  specified.  To  separate  the  different  domains,  we 
introduce  a  level  set  function  (j)  defined  as: 

!(f>  <  0  for  x  €  D~, 

4>  >  0  for  x  €  D+, 

4>  =  0  for  x  6  T. 

A  convenient  choice  that  ensures  numerical  robustness  is  to  define  4>  as  the 
signed  distance  function  to  the  interface.  The  level  set  function  is  also  used 
to  identify  the  location  of  the  interface  to  high  order  accuracy  as  will  be 
discussed  throughout  this  paper. 

The  discretization  of  the  Laplace  operator,  including  the  special  treat¬ 
ments  needed  at  the  interface,  is  performed  in  a  dimension  by  dimension 
fashion.  Therefore,  without  loss  of  generality,  we  first  describe  the  dis¬ 
cretization  in  one  spatial  dimension. 
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2.1  ID  Laplace  Equation 

Consider  the  Laplace  equation  in  one  spatial  dimension,  i.e.  Txx  =  /.  The 
computational  domain  is  discretized  into  cells  of  size  Ax  with  the  grid  nodes 
Xi  located  at  their  centers.  The  solution  to  the  Laplace  equation  is  computed 
at  the  grid  nodes  and  is  written  as  Tt  =  T(xj).  We  consider  the  standard 
fourth  order  discretization: 


(T'xx)i 


1  rT.  |  4  rji  5t  i  4t 

2  '  3 1  2^^  '  3-M+l 

Ax2 


12  Ti+ 2 


(2) 


For  each  unknown,  T*,  equation  (2)  is  used  to  fill  in  one  row  of  a  matrix 
creating  a  linear  system  of  equations.  This  discretization  is  valid  if  all  the 
nodal  values  used  belong  to  the  same  domain,  but  needs  to  be  modified 
otherwise.  For  example,  suppose  the  interface  location  xj  is  located  in 
between  the  nodes  Xi  and  Xj+i  (see  figure  1),  and  suppose  that  we  seek  to 
write  the  equation  satisfied  by  T*.  Since  the  solution  is  not  defined  across  the 
interface,  we  need  valid  values  for  \  and  T,;+ 2  that  emulate  the  behavior  of 
the  solution  defined  to  the  left  of  the  interface.  We  achieve  this  by  defining 
‘ghost  values’  T^1  and  T^_2  constructed  by  extrapolating  the  values  of  T 
across  the  interface.  The  discretization  is  then  rewritten  as 


( Txx)i 


_ Lt1  _i_  ^rr  _ 

l2-t*-2  “I"  1  9-tl 


Ax2 


4  TG 

3L.+1 


It iG 
121i+2 


(3) 


More  precisely,  we  first  construct  an  interpolant  T{x)  of  T(x )  on  the 
left  of  the  interface  with  T(0)  =  Tj,  and  then  we  define  =  T(Ax)  and 
TB_ 2  =  T{2Ax).  Figure  1  illustrates  the  definition  of  the  ghost  cells  in  the 
case  of  the  linear  extrapolation.  In  this  paper  we  consider  constant,  linear, 
quadratic  and  cubic  extrapolation  defined  by: 

Constant  extrapolation:  Take  T(x)  =  d  with 

•  d  =  2>. 

Linear  extrapolation:  Take  T(x)  =  cx  +  d  with 

•  T(  0)  =  Tu 

•  T{6Ax)  =  Tj. 

Quadratic  extrapolation:  Take  T(x)  =  bx2  +  cx  +  d  with 

•  T(-Ax)  =  Tj_  1, 
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•  T(  0)  =  Ti, 

•  T{9Ax)  =  Tj. 

Cubic  extrapolation:  Take  Z(x)  =  ax 3  +  bx 2  +  cx  +  d  with 

•  Z(— 2Ax)  =  Ti_  2, 

•  Z(-Ax)  =  i, 

•  f(0)  =  Ti , 

•  Z(0Ax)  =  Tj. 

In  these  equations  0  =  (xj  —  Xi)/Ax. 

Similarly,  if  we  were  solving  for  the  domain  ,  the  equation  satisfied  by 
Tj+i  requires  the  definition  of  the  ghost  cells  Tf*  and  T-_j .  In  this  case,  we 
write  Zp  =  Z(Ax)  and  Zicl1  =  Z(2Ax)  with  the  definition  for  T  modified 
as  follows:  8  is  replaced  by  1  —  9,  Ti  is  replaced  by  Tl+\.  Tt-\  is  replaced  by 
Ti+ 2  and  Z;_2  is  replaced  by  Ti+3. 

The  construction  of  Z  is  both  limited  by  the  number  of  points  in  the  do¬ 
main  and  by  how  close  the  interface  is  to  a  grid  node.  The  latter  restriction 
stems  from  the  fact  that  as  9  — >  0,  the  behavior  of  the  interpolant  deterio¬ 
rates.  We  found  that  a  good  rule  of  thumb  is  to  shift  the  interpolation  to  be 
centered  one  grid  point  to  the  left  when  9  <  Ax,  e.g.  in  the  case  of  a  linear 
extrapolation,  we  use  the  conditions  Z(0)  =  Z*_i  and  T( Ax  +  9 Ax)  =  Tj 
instead  of  Z(0)  =  Ti  and  T{6 Ax)  =  Tj.  Then  the  ghost  nodes  are  defined 
as  Zj®1  =  Z( 2 Ax)  and  Z^2  =  Z( 3 Ax).  Higher  order  extrapolation  follows 
similarly.  Finally,  we  note  that  one  can  lower  the  degree  of  the  interpolant 
in  order  to  preserve  the  pentadiagonal  structure  of  the  linear  system  in  the 
case  where  the  stencil  has  been  shifted. 

For  the  constant  and  the  linear  extrapolations,  the  matrix  entries  to  the 
right  of  the  diagonal  for  the  z-th  row  and  to  the  left  of  the  diagonal  for  the 
(i  +  l)-th  row  are  equal  to  zero,  yielding  a  symmetric  linear  system.  This  al¬ 
lows  for  the  use  of  fast  iterative  solvers  such  as  the  preconditioned  conjugate 
gradient  (PCG)  method  (see,  e.g.  [28]).  Moreover  the  corresponding  matrix 
is  strictly  diagonally  dominant,  and  therefore  nonsingular.  For  higher  or¬ 
der  extrapolations,  the  linear  system  is  nonsymmetric  and  not  necessarily 
strictly  diagonally  dominant,  but  we  can  still  develop  high  order  accurate 
methods.  The  overall  accuracy  for  Z  and  the  nature  of  the  resulting  linear 
system  is  determined  by  the  degree  of  the  interpolation  which  is  summarized 
in  table  1.  Finally,  we  note  that  [21]  designed  a  method  that  also  yields  a 
nonsymmetric  linear  system,  but  which  is  only  second  order  accurate. 
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Degree  of  Extrapolation 

Order  of  Accuracy 

Linear  System 

Constant 

First 

Symmetric 

Linear 

Second 

Symmetric 

Quadratic 

Third 

Nonsymmetric 

Cubic 

Fourth 

Nonsymmetric 

Table  1:  Order  of  accuracy  and  nature  of  the  linear  system. 


We  illustrate  our  methodology  with  the  following  example.  Let  17  =  [0, 1] 
with  an  exact  solution  of  T  =  x 5  —  x3  +  12.x2  —  2.5x  +  2  on  fi-.  We  define 
(j)  =  x  —  .5  (so  that  the  interface  never  falls  on  a  grid  node).  Dirichlet 
boundary  conditions  are  enforced  on  the  <917  using  the  exact  solution.  Tables 
2-5  give  the  error  between  the  numerical  solution  and  the  exact  solution  in 
the  L 1  and  L°°-norms.  These  same  results  are  also  presented  on  a  log-log 
plot  in  figure  2,  where  the  open  symbols  represent  the  error  in  the  L°°-norm 
and  the  solid  lines  represent  the  least  square  fit.  These  results  illustrate  the 
first,  second,  third  and  fourth  order  accuracy  in  the  case  of  constant,  linear, 
quadratic  and  cubic  extrapolation,  respectively.  We  use  a  PCG  method  with 
an  incomplete  Cholesky  preconditioner  for  the  symmetric  linear  systems,  and 
the  BiCGSTAB  method  (see  e.g.  [28])  otherwise. 


Number  of  Points 

L1  —  error 

order 

L°°  —  error 

order 

16 

1.307  x  10-1 

— 

2.369  x  10"1 

— 

32 

6.248  x  10"2 

1.06 

1.196  x  10"1 

.98 

64 

3.057  x  10-2 

1.03 

6.018  x  10-2 

.99 

128 

1.512  x  10~2 

1.02 

3.020  x  10"2 

1.00 

Table  2:  Accuracy  results  for  the  one  dimensional  Laplace  equation  with 
ghost  cells  defined  by  constant  extrapolation. 


2.2  2D  Laplace  Equation 


The  methodology  discussed  in  section  2.1  extends  naturally  to  two  and  three 
spatial  dimensions.  For  example,  in  the  case  of  two  spatial  dimensions,  we 
solve 


TXx  +  Tyy  —  f  ■ 
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Number  of  Points 

L1  —  error 

order 

L°°  —  error 

order 

16 

4.456  x  1CT3 

— 

8.463  x  10~3 

— 

32 

1.013  x  10~3 

2.13 

2.045  x  10“3 

2.05 

64 

2.417  x  10“4 

2.06 

5.031  x  10~4 

2.02 

128 

5.901  x  10“5 

2.03 

1.247  x  10“4 

2.01 

Table  3:  Accuracy  results  for  the  one  dimensional  Laplace  equation  with 
ghost  cells  defined  by  linear  extrapolation. 


Number  of  Points 

L1  —  error 

order 

L°°  —  error 

order 

16 

2.168  x  10”5 

— 

5.197  x  10"5 

— 

32 

3.084  x  10"6 

2.81 

7.532  x  10~6 

2.78 

64 

4.013  x  10~7 

2.94 

9.971  x  10“7 

2.91 

128 

5.095  x  10~8 

2.98 

1.278  x  10“7 

2.96 

Table  4:  Accuracy  results  for  the  one  dimensional  Laplace  equation  with 
ghost  cells  defined  by  quadratic  extrapolation. 
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and  for  cells  cut  by  the  interface,  ghost  values  are  defined  by  extrapolating 
the  value  of  T  across  the  interface  as  described  in  section  2.1.  In  two  spatial 
dimensions,  the  definition  of  the  ghost  cells  involves  9X  and  9y,  i.e.  the 
cell  fractions  in  the  x  and  y  direction,  respectively.  These  quantities  are 
evaluated  as  follows.  Consider  a  grid  node  ( Xi,yj )  in  the  neighborhood  of 
the  interface.  We  first  construct  a  cubic  interpolant  c j>x  of  4>  in  the  x-direction 
and  find  the  interface  location  xj  by  solving  (f>x(xi)  =  0.  Then  we  define 
9X  =  | Xi  —  xj  |/Ax.  The  procedure  to  find  9y  is  similar.  We  emphasize  that 
the  numerical  discretization  of  Txx  is  independent  from  that  of  Tyy ,  making 
the  procedure  trivial  to  extend  to  two  and  three  spatial  dimensions. 

We  illustrate  the  order  of  accuracy  with  the  following  example.  Let  12  = 
[— 1, 1]  x  [—1, 1]  with  an  exact  solution  of  T  =  sin(7rx)  +  sin(7ry)  +cos(7rx)  + 


Number  of  Points 

L1  —  error 

order 

L°°  —  error 

order 

16 

1.502  x  10_(i 

— 

8.519  x  10~6 

— 

32 

8.416  x  10"8 

4.15 

5.401  x  10"7 

3.97 

64 

4.867  x  10“9 

4.11 

3.378  x  10“8 

3.99 

128 

2.936  x  10~10 

4.05 

2.109  x  10~9 

4.00 

Table  5:  Accuracy  results  for  the  one  dimensional  Laplace  equation  with 
ghost  cells  defined  by  cubic  extrapolation. 

cos(7ry)  +  x6  +  ye  in  H-.  The  interface  is  parameterized  by  (x(a),y(a)) 
where: 

I  x(a)  =  .02^/5  +  (.5  +  .2  sin(5a))  cos(a), 

\  y(a )  =  .02v/5  +  (.5  +  .2  sin(5a))  sin(a), 

with  a  €  [0,  2tx].  Figure  3  depicts  the  solution  on  a  64  x  64  grid,  and  figure  4 
illustrates  the  accuracy  in  the  L°°-norm.  Note  that  on  irregular  domains,  the 
number  of  available  grid  nodes  within  the  domain  might  limit  the  extrapola¬ 
tion  to  a  lower  degree  for  some  grid  resolutions.  This  partially  explains  the 
‘oscillatory’  nature  of  the  accuracy  results  in  the  graph  of  figure  3.  However, 
the  slopes  of  the  least  square  fits  are  still  in  accordance  with  first,  second, 
third  and  fourth  order  accuracy  for  the  constant,  linear,  quadratic  and  cubic 
extrapolation. 

3  Heat  Equation 

Consider  a  Cartesian  computational  domain,  H  C  Rn,  with  exterior  bound¬ 
ary,  90,  and  a  lower  dimensional  interface,  T,  that  divides  the  computational 
domain  into  disjoint  pieces,  and  0+.  The  Heat  equation  is  written  as 

Tt(x )  =  A T(x),  x  €  0~,  (4) 

where  T{x)  is  assumed  to  be  smooth  on  0~.  Dirichlet  boundary  conditions 
are  specified  on  T. 

Explicit  time  discretization  schemes  are  impractical  in  the  case  of  arbi¬ 
trary  domains,  because  they  suffer  from  stringent  time  step  restrictions.  For 
example  in  one  spatial  dimension,  we  must  impose  a  time  step  restriction  of 
O  (92Ax2)  with  0  <  9  =  (xj— Xj) /  Ax  <  1  for  cells  cut  by  the  interface.  Since 
9  can  be  arbitrarily  small,  explicit  schemes  are  prohibitively  computation¬ 
ally  expensive.  Although  one  could  remesh  the  domain  to  keep  9  reasonable, 


9 


in  the  case  of  a  moving  interface  this  would  require  remeshing  every  time 
the  value  of  0  gets  below  an  acceptable  threshold.  Thus,  we  use  implicit 
time  discretization.  In  particular,  we  choose  the  Crank-Nicholson  scheme 
and  impose  a  time  step  restriction  of  At  =  cAx2  with  0  <  c  <  1  in  order  to 
obtain  a  fourth  order  accurate  discretization  in  time.  However,  we  note  that 
Backward  Differentiation  Formulae  or  implicit  Runge-Kutta  schemes  could 
be  use  instead  in  order  to  relax  the  time  step  restriction  to  At  =  cAx.  For 
more  on  numerical  methods  for  ordinary  differential  equations,  see  [20]. 

The  Crank-Nicholson  scheme  can  be  written  as 

(l+^An{ A))  Tn, 

where  An(  A)  and  Hn+1(  A)  represent  the  spatial  approximation  of  the  Laplace 
operator  at  time  tn  and  tn+1,  respectively.  The  spatial  discretization  is  per¬ 
formed  in  a  dimension  by  dimension  fashion  and  resembles  that  of  section  2. 
More  precisely,  we  first  evaluate  the  right  hand  side  fn  =  (V  +  Tn 

at  time  tn  using  the  methodology  of  section  2  to  define  the  ghost  cells  in  a 
dimension  by  dimension  fashion.  Then,  we  solve 

(i  -  ^~An+1{ A)j  Tn+1  =  fn. 

The  Laplace  operator  at  time  tn+1  is  discretized  along  the  lines  of  section  2 
as  well.  Each  equation  is  used  to  fill  one  row  of  a  linear  system  that  is  then 
solved  with  an  iterative  solver. 

Since  the  discretization  is  performed  in  a  dimension  by  dimension  fash¬ 
ion,  we  first  present  the  one  dimensional  case. 

3.1  ID  Heat  Equation 

In  one  spatial  dimension,  we  discretize  the  heat  equation  as 

Tn+ 1  -  ^ An+1(TXX )  =  Tn  +  ^An(Txx), 

where  An(Txx)  and  An+1(TXX )  are  the  fourth  order  approximations  of  Txx 
at  time  tn  and  tn+1,  respectively.  The  discretization  of  the  heat  equation  is 
performed  in  two  steps  and  depends  heavily  on  that  of  the  Laplace  operator 
described  in  section  2.1.  First,  we  approximate  Txx  with  the  fourth  order 
accurate  discretization  of 

_ )^rpn  4 nnn  bn~<n  4nrn  1  rpn 

(rpn  \  ^  12 1  i—2  '  3 z— 1  2 1  i  f  3^+1  121  i+2 

K  xxh  ~  Arc2 


(l-  ^An+\A)^j  Tn+1  = 


10 


The  special  treatment  needed  for  grid  nodes  neighboring  the  interface  is 
performed  as  described  in  section  2.1,  using  the  values  at  the  interface  at 
time  tn.  We  then  evaluate  f  n  =  Tn  +  ^  An{Txx )  and  are  left  to  solve 


Tn+l  _  _An+l{Txx)  =  fn  (5) 

We  again  employ  the  standard  fourth  order  discretization  to  approximate 
Txx  at  time  tn+1 


1  rrin-\-l  I  4rpn-\-l  _  5_rpn-\- 1  Arpn-\-l  _  1  rrm+l 

12 1  i—2  '  S1i—1  2 '  3^+1  12^+2 

Ax2 


(6) 


For  each  unknown,  T”+1,  equations  (5)  and  (6)  are  used  to  fill  in  one  row 
of  a  matrix  creating  a  linear  system  of  equations.  The  treatment  of  the  grid 
nodes  in  the  neighborhood  of  the  interface  is  again  based  on  defining  ghost 
node  values  and  uses  the  values  of  the  temperature  at  the  interface  at  time 
tn+1. 

For  the  remainder  of  this  paper,  we  focus  on  designing  a  fourth  order 
accurate  method  for  the  heat  equation  and  a  third  order  accurate  method 
for  the  Stefan  problem.  Therefore,  we  present  only  the  results  for  these 
accuracy  tests,  although  we  have  checked  that  one  obtains  first,  second, 
third  and  fourth  order  accuracy  for  the  constant,  linear,  quadratic  and  cubic 
extrapolation,  respectively.  The  nature  of  the  linear  system  and  the  order 
of  accuracy  is  the  same  as  that  of  the  Laplace  operator  (see  Table  1). 

Consider  the  following  example.  Let  fl  =  [—1,1]  with  an  exact  solution 
of  T  =  e-71"2*  cos(7rx)  on  f 1~ .  We  take  (j)  =  x  —  .313  (thus  0  <  9  <  1).  Dirich- 
let  boundary  conditions  are  enforced  on  the  dVL  using  the  exact  solution, 
and  the  final  time  is  t  =  \  ji r2.  The  ghost  values  are  defined  with  a  cubic 
extrapolation  of  T  across  the  interface.  Figure  5  illustrates  the  fourth  order 
accuracy  in  the  L°°-norm. 


3.2  2D  Heat  Equation 

The  algorithm  described  above  extends  readily  to  two  and  three  spatial  di¬ 
mensions.  The  approximations  of  Txx  and  Tyy  are  performed  independently 
making  the  procedure  trivial  to  implement.  Consider  the  following  example. 
Let  Q  =  [—1,1]  x  [—1,1]  with  an  exact  solution  of  T  =  e~2t  sin(x)  sin(y)  in 
D-.  The  interface  is  parameterized  by  (x(a),y(a))  where: 

f  x(a)  =  .02v/5  +  (.5  +  .2  sin(5a))  cos(a), 

1  y(a)  =  .02v/5  +  (.5+.2  sin(5a))  sin(a), 
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with  a  e  [0,  2ir\.  Figure  6  depicts  the  solution  on  a  64  x  64  grid  and  figure 
7  demonstrates  the  fourth  order  accuracy  in  the  L°°-norm. 

4  Stefan  Problem 

Consider  again  a  Cartesian  computational  domain,  O  C  Rn,  with  exterior 
boundary,  <90,  and  a  lower  dimensional  interface,  T,  that  divides  the  com¬ 
putational  domain  into  disjoint  pieces,  and  0+.  The  Stefan  problem  is 
written  as 

(  Tt(x)  =  DAT(x),  fed, 

<  T(x)  =  0,  fer, 

{  Vn  =  ~  [DVT]  r  •  n, 

where  D  is  the  diffusion  coefficient,  assumed  in  this  work  to  be  constant  on 
each  subdomain  (but  possibly  discontinuous  across  the  interface).  T( x)  is 
assumed  to  be  smooth  on  each  disjoint  subdomain,  O-  and  0+,  but  may 
have  a  kink  at  the  interface  F.  Diri chief  boundary  conditions  are  specified 
on  <9D. 

We  need  to  both  discretize  the  heat  equation  and  evaluate  the  velocity 
at  the  interface.  The  added  complexity  for  the  Stefan  problem  stems  from 
the  fact  that  the  interface  is  evolving  in  time.  We  keep  track  of  the  interface 
evolving  under  the  velocity  field  V  =  (it,  v,  w)  by  solving  the  advection 
equation 


4>t  +  V  ■  V(f>  =  0. 

The  velocity  components  are  defined  by  the  x,  y  and  z  projections  of  the 
jump  in  the  temperature  gradient,  i.e.  ( u,v,w )  =  {[DTX\ r,  [DTy] p,  [DTZ] p). 
The  level  set  advection  equation  is  discretized  with  a  HJ-WENO  scheme 
[12],  see  also  [22,  13].  For  more  details  on  the  level  set  method,  see  e.g. 
[24,  30]. 

Consider  the  Crank-Nicholson  framework: 

Tn+i  _  ^An+1(A)Tn+1  =  Tn  +  ^-An(A)Tn,  (7) 

and  let  the  temperature  be  defined  at  time  tn.  The  algorithm  to  solve  the 
Stefan  problem  is: 

1.  Extrapolate  Tn  in  the  normal  direction. 

2.  Discretize  fn  =  Tn  +  ^An(A)Tn. 
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3.  Evolve  the  level  set  function  for  one  time  step  At. 

4.  Assemble  and  solve  the  linear  system  for  Tn+1. 

5.  Repeat  1-4  until  done. 


4.1  Extrapolation  in  the  Normal  Direction 


We  follow  along  the  lines  of  section  3  discretizing  the  heat  equation  with  the 
Crank-Nicholson  scheme  by  writing  the  discretization  of  the  Laplace  opera¬ 
tor  at  time  tn  and  evaluating  the  right  hand  side  fn  =  Tn  +  ^ An(A)Tn . 
However,  the  moving  domain  provides  added  difficulty.  Figure  8  depicts 
a  two  dimensional  example.  As  the  interface  moves  from  its  position  at 
time  tn  to  its  new  location  at  time  tn+1,  new  grid  nodes  are  added  to  IE 
(depicted  for  example  by  the  dark  node  in  the  figure).  Since  we  need  to 
evaluate  the  right  hand  side  at  these  new  nodes,  valid  values  for  Tn  must  be 
defined  there.  Moreover,  since  the  interface  evolves  in  the  normal  direction, 
valid  values  of  Tn  are  needed  in  more  than  just  the  Cartesian  directions. 
Therefore,  a  high  order  extrapolant  must  be  defined  in  the  normal  direction 
at  the  interface  and  such  a  procedure  must  be  straightforward  to  implement 
in  one,  two  and  three  spatial  dimensions. 

High  order  extrapolation  in  the  normal  direction  is  performed  in  a  series 
of  steps,  as  proposed  in  [1] .  For  example,  suppose  that  we  seek  to  extrapolate 
T  from  the  region  where  <f>  <  0  to  the  region  where  (f>  >  0.  In  the  case  of 
cubic  extrapolation,  we  first  compute  Tnnn  =  V  (VT  ■  nj  ■  nj  ■  ft  in  the 
region  cf>  <  0  and  extrapolate  it  across  the  interface  in  a  constant  fashion  by 
solving  the  following  partial  differential  equation 


+  H(4>  +  band)VTnnn  ■  n  =  0, 

OT 

where  H  is  the  Heaviside  function  and  band  accounts  for  the  fact  that  Tnnn 
is  not  defined  in  the  region  where  (f)  >  —band.  For  example,  we  set  band  = 
2 \J dx2  +  dy2  in  the  case  where  Tnnn  is  computed  by  central  differencing. 
Then,  the  value  of  T  across  the  interface  is  found  by  solving  the  following 
three  partial  differential  equations.  First,  solve 


dTn 


dr 


+ 


H(<j>)  (VT, 


—  T 

nn  nnn 


=  o, 


defining  Tnn  in  such  a  way  that  its  normal  derivative  is  equal  to  Tnnn.  Then 
solve 


dTn 

dr 


+  H(4>)  (VTn  —  Tn„J  —  0, 
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defining  Tn  in  such  a  way  that  its  normal  derivative  is  equal  to  Tnn .  Finally 
solve 

—  +  H(</>)  (VT  -  Tn)  =  0, 

defining  T  in  such  a  way  that  its  normal  derivative  is  equal  to  Tn.  These 
equations  are  solved  in  fictitious  time  r  for  a  few  iterations  (typically  15), 
since  we  only  seek  to  extrapolate  the  values  of  T  in  a  narrow  band  of  a  few 
grid  cells  around  the  interface. 

Figure  9  illustrates  cubic  extrapolation.  This  example  is  taken  from  [1]. 
Consider  a  computational  domain  17  =  [— it,  7r]  x  [— 7t,  7t]  separated  into  two 
regions:  17“  defined  as  the  interior  of  a  disk  with  center  at  the  origin  and 
radius  two,  and  its  complementary  17+.  The  function  T  to  be  extrapolated 
from  17~  to  17+  is  defined  as  T  =  cos(x)  sin(y)  for  x  G  !7-.  Figure  9  (top) 
illustrates  contours  of  T  after  it  has  been  extrapolated  across  the  interface, 
and  figure  9  (bottom)  depicts  contours  of  the  exact  solution  for  comparison. 
For  the  sake  of  clarity,  we  have  extrapolated  T  in  the  entire  region  in  this 
example,  but  in  practice  the  extrapolation  is  performed  only  in  a  neigh¬ 
borhood  of  the  interface.  The  extrapolation  is  fourth  order  accurate  in  the 
L^-norm  as  demonstrated  in  figure  10. 

At  the  end  of  the  high  order  extrapolation  procedure  we  have  Tn  at  all 
grid  nodes  near  the  interface,  so  that  if  the  interface  moves  across  new  grid 
nodes,  we  can  still  evaluate  Tn  +  ^An(A)Tn  and  proceed  as  in  section  3 
to  assemble  the  linear  system. 

4.2  Discretization  of  the  Velocity 

The  method  described  above  can  be  used  to  obtain  a  fourth  order  accurate 
temperature.  Consequently,  the  velocity  will  only  be  third  order  accurate 
as  it  involves  gradients  of  the  temperature  .  Not  surprisingly,  this  limits  the 
overall  order  of  accuracy  to  three.  For  example,  let  17  =  [0, 1]  with  an  exact 
solution  of  T  =  e — i  Qn  17”  with  4>  =  x— .5  at  t  =  0.  Dirichlet  boundary 
conditions  are  enforced  on  <917  using  the  exact  solution.  We  compute  the 
solution  at  time  ffinai  =  .25.  We  use  the  the  Crank-Nicholson  scheme  with 
A t  ~  Ax2  for  fourth  order  time  accuracy,  and  use  cubic  extrapolation  to 
define  the  ghost  values.  In  a  sequence  of  tests,  we  use  the  exact  interface 
velocity  perturbed  by  an  0(Ax3),  0(Ax2)  and  O(Ax)  amount.  Figure  11 
illustrates  the  fact  that  a  third,  second  and  first  order  accurate  velocity 
produces  overall  results  of  the  same  order. 

The  examples  above  demonstrate  that  a  O(Ax)  perturbation  in  the  ve¬ 
locity  forces  the  solution  to  be  first  order  accurate  in  the  case  of  the  Stefan 
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problem,  regardless  of  the  order  of  accuracy  of  the  discretization  of  the  heat 
equation  operator.  In  [9],  the  velocity  could  be  at  most  first  order  accurate, 
since  it  was  computed  as  a  derivative  of  a  second  order  accurate  solution. 
As  a  consequence,  the  authors  had  the  leeway  to  compute  the  jumps  in 
Tx  (respectively  Ty  and  Tz)  at  the  point  where  the  interface  crosses  the  x 
(respectively  y  and  z)  axis.  That  short  cut  in  the  velocity  computation  pro¬ 
duces  a  simple  discretization,  but  it  cannot  readily  be  extended  to  higher 
order  accuracy.  Moreover,  since  the  level  set  moves  in  the  normal  direction, 
the  velocity  should  be  constant  in  the  normal  direction.  That  is,  the  velocity 
at  a  grid  node  near  the  interface  should  be  equal  to  that  of  the  closest  point 
on  the  interface.  In  addition,  in  order  to  construct  an  overall  third  order 
accurate  scheme,  we  require  a  third  order  accurate  velocity. 

Suppose  that  we  construct  a  cubic  interpolation  T  of  the  temperature 
around  the  interface,  and  that  the  interface  position  xj  is  known  to  third 
order  accuracy.  Then  the  velocity  is  defined  as  ([Tx]r,  [T^Jr,  [Tz]r).  The 
construction  of  T  is  straightforward  once  the  temperature  values  have  been 
extrapolated  across  the  interface  as  described  in  section  4.1,  since  we  then 
have  valid  values  for  the  temperature  of  each  domain  in  both  (f>  >  0  and 
(f)  <0. 

4.2.1  Finding  the  Closest  Point 

There  are  many  ways  of  finding  the  closest  point  to  an  implicitly  defined 
interface.  Here,  we  follow  the  work  of  [3]  since  it  is  based  on  bi-cubic  interpo¬ 
lation  that  fits  well  into  our  framework.  Given  the  level  set  function,  one  can 
identify  the  cells  crossed  by  the  interface  by  simply  checking  the  sign  change 
of  (j>.  For  each  such  cell  C  with  vertices  ( Xi,yj ),  (xi+i,yj),  (xi,yj+ 1)  and 
(xj+i,  yj+i),  we  construct  a  cubic  interpolation  (f>  of  <j>  using  the  grid  nodes 
of  the  3x3  cells  centered  at  C.  For  each  grid  node  P  =  ( x^yj )  near  the 
interface,  we  seek  to  find  the  closest  point  on  the  set  S  =  |x  €  £l\(f)(x)  =  0  j. 
[3]  notes  that  such  a  point  x}  must  satisfy  the  following  two  equations: 

cj)(xi)  =  0,  (8) 

V<f>  x  (P  —  xj)  =  0,  (9) 

accounting  for  the  fact  that  xj  must  be  on  S  and  that  the  normal  at  this 
point  must  be  aligned  with  xj  —  P.  Then  [3]  proposes  an  iterative  scheme 
starting  with  xp]  =  P  and  solving  simultaneously  equations  (8)  and  (9)  with 
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a  Newton-type  algorithm: 
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Convergence  is  assumed  when  \/||<5i||2  +  ||d2||2  <  10-3  Ax  Ay.  Typically,  five 
iterations  are  sufficient  to  find  the  closest  point  to  third  order  accuracy  in 
the  case  where  the  interface  is  smooth.  However,  we  use  20  iterations  just 
to  be  safe.  Although  convergence  is  not  guaranteed,  we  did  not  encounter 
any  problems  in  our  computations.  The  reader  is  referred  to  [3]  for  more 
details. 


4.2.2  A  Note  on  the  Reinitialization  Equation 

For  the  sake  of  robustness  in  the  numerics,  it  is  important  to  keep  the  values 
of  (j)  close  to  those  of  a  signed  distance  function,  i.e.  |V<£|  =  1.  Since  the 
Fast  Marching  Method  [34,  29]  is  only  first  order  accurate,  and  an  0( Ax) 
perturbation  of  the  interface  leads  to  a  first  order  accurate  temperature, 
we  cannot  use  this  method.  Another  way  of  reinitializing  <f>  is  to  solve  the 
reinitialization  equation  introduced  in  [33] 

<f>T  +  S((t>o)  (|v</>|  -  1)  =  0  (10) 

for  a  few  steps  in  fictitious  time,  r.  This  equation  is  traditionally  discretized 
with  the  HJ-WENO  scheme  of  [12]  in  space,  because  it  yields  less  noisy 
results  when  computing  quantities  such  as  the  curvature.  However,  we  note 
that  this  method  is  only  second  order  accurate  as  depicted  in  figure  12.  This 
is  due  to  the  fact  that  the  reinitialization  equation  is  a  Hamilton- Jacobi  type 
equation  with  discontinuous  characteristics  emanating  from  the  interface. 
Since  a  0(  Ax2)  perturbation  of  the  interface  leads  to  a  second  order  accurate 
solution,  we  do  not  use  equation  (10)  to  reinitialize  <j>.  Instead,  the  procedure 
detailed  in  section  4.2.1  gives  the  distance  to  the  interface  to  third  order 
accuracy  for  grid  nodes  in  the  neighborhood  of  the  interface. 

4.2.3  A  Simple  Example  in  One  Spatial  Dimension 

Let  12  =  [0, 1]  and  (j)  =  x  —  .5  at  t  =  0.  We  consider  an  exact  solution 
of  T  =  et_a:+'5  —  1  on  IW.  Dirichlet  boundary  conditions  are  enforced 
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on  the  dQ  using  the  exact  solution,  and  we  compute  the  solution  at  time 
t Anal  =  -25.  We  use  the  Crank-Nicholson  scheme  with  At  ~  Ax2,  and  cubic 
extrapolation  to  define  the  ghost  cells.  Figure  13  demonstrates  the  fourth 
order  accuracy  obtained  when  using  the  exact  interface  velocity  and  the 
third  order  accuracy  obtained  when  using  the  computed  interface  velocity. 

4.3  Time  Discretization 

In  the  example  in  section  4.2.3,  the  interface  velocity  is  constant  in  time, 
i.e.  [Tx]  r  =  1.  Therefore,  this  example  focused  on  the  spatial  discretization 
and  velocity  computation,  but  was  not  sensitive  to  the  time  discretization 
details. 

Consider  the  Frank  spheres  solution  in  one  spatial  dimension.  Let  14  = 
[—1,1]  with  Dirichlet  boundary  conditions  at  the  domain  boundaries.  The 
Frank  spheres  solution  in  one  spatial  dimension  describes  a  slab  of  radius 
R(t)  =  So\/i  parameterized  by  Sq.  The  exact  solution  takes  the  form 

JO  s  <  S0, 

~  {  r~  (i  -  M)  8  >  s»- 

where  s  =  \x\/\ft.  T, 'x  and  Sq  are  related  by  the  jump  condition  Vn  = 
— D  [VT]r  •  n.  In  one  spatial  dimension  F(s)  =  erfc(s/2),  with  erfc(z)  = 
2  /2°°  dtfypK.  We  choose  tinitiai  =  1  and  =  —.5  obtaining  an  initial 
radius  of  Sq  «  .86.  Initially,  we  set  <j>  =  \x\  —  Sq  and  compute  the  solution 
at  time  tfinai  =  1.5  with  the  Crank-Nicholson  scheme.  Since  the  velocity  is 
third  order  accurate  (hence  we  expect  a  third  order  accurate  method),  we 
choose  a  time  step  restriction  of  At  ~  Ax3/2  to  obtain  a  third  order  accurate 
scheme  in  time  as  well.  However,  the  time  discretization  presented  so  far 
yields  results  that  are  only  second  order  accurate  for  time  varying  velocities 
as  demonstrated  in  figure  14. 

This  lower  accuracy  stems  from  the  lack  of  consistency  in  the  definition 
of  Vn+1.  For  example,  consider  approximating 

f  =  n*). 

with  the  Crank-Nicholson  scheme.  To  evolve  from  time  tn  to  time  tn+1, 
we  perform  the  following  three  steps: 

1.  Use  Vn((j)n)  to  evolve  ^>n  to  (frtemp  an  Euler  step. 

2.  Use  Vn+1  {(Pt^rip)  evolve  to  <f)n+ 2  with  an  Euler  step. 
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3.  Define  (j)n+1  =  (<pn  +  (j)n+2)/ 2. 

In  the  case  of  the  Stefan  problem,  the  velocity  at  time  tn+1  is  defined  from 
the  jump  in  the  temperature  gradient  at  the  interface  at  time  tn+1,  and 
thus  needs  to  be  consistent.  That  is,  if  Vn+l  =  Vn+1  ((j)n+1) ,  then  the 
yn+\  from  g^ep  2  above  needs  to  be  consistent  with  the  (j)n+1  computed  in 
step  3.  To  solve  this  problem  we  iterate  steps  2  and  3  in  order  to  guaran¬ 
tee  that  the  velocity  at  time  tn+ 1  is  consistent,  i.e  Vn+l  =  Vn+1((f>n+1)  = 
V  ((</>n  +  4>n+2)  /2).  We  iterate  until  the  magnitude  of  the  error  in  this  last 
equation  is  less  than  10“8  noting  that  very  few  iterations  are  needed.  Figure 
15  demonstrates  that  this  time  discretization  yields  a  third  order  accurate 
solution. 

4.4  Example  in  Two  Spatial  Dimensions 

Consider  the  Stefan  problem  in  a  domain  [—1,1]  x  [—1,1]  with  Dirichlet 
boundary  conditions  at  the  domain  boundaries.  The  Frank  spheres  in  two 
spatial  dimensions  describes  a  disk  of  radius  R(t)  =  SoVt  parameterized  by 
Sq.  The  exact  solution  takes  the  form 

f  °  s  <  50, 

\  roo  (l  -  Fpo))  s  >  So, 

where  s  =  \x\/y/t.  T, 'x  and  So  are  related  by  the  jump  condition  Vn  = 
— D  [VT]r  •  n.  In  two  spatial  dimensions,  F{s)  =  Ei(s2/4)  with  Ei(z)  = 
/2°°  e_t /tdt.  We  take  tinitiai  =  1  and  So  =  .5  ((f)  =  |x|  —  So).  This  defines 
Too  ~  —.15.  We  compute  the  solution  at  time  ffinai  =  2.89  with  the  ghost 
cells  defined  via  cubic  extrapolation.  The  Frank  spheres  problem  being 
ill-posed,  we  choose  the  final  time  large  enough  to  allow  the  interface  to 
cross  a  large  amount  of  grid  cells  (about  50),  demonstrating  the  built-in 
regularization  inherent  to  level  set  methods.  Table  6  presents  the  accuracy 
results. 

For  the  sake  of  comparison,  table  7  offers  the  convergence  results  ob¬ 
tained  with  the  symmetric  discretization  from  [9] .  We  note  that  the  present 
method  and  16  grid  nodes  yields  the  same  accuracy  as  in  [9]  with  128  grid 
nodes.  Likewise,  [9]  would  require  1080  points  to  obtain  the  same  accuracy 
as  the  present  method  with  32  points.  This  may  have  a  significant  impact, 
especially  in  three  spatial  dimensions.  In  fact,  even  when  utilizing  adap¬ 
tive  grid  refinement,  this  newly  proposed  method  can  drastically  lower  the 
number  of  grid  nodes  needed  to  represent  thin  dendrites  while  retaining  the 
desired  accuracy.  Figure  18  illustrates  the  comparison  between  the  method 
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Number  of  Points 

L1  —  error 

order 

L°°  —  error 

order 

16  x  16 

3.204  x  10-4 

— 

1.032  x  10~3 

— 

32  x  32 

1.092  x  10~5 

4.87 

6.954  x  10~5 

3.89 

64  x  64 

7.369  x  10“7 

3.89 

3.482  x  10~6 

4.32 

128  x  128 

8.836  x  10“8 

3.06 

3.149  x  10~7 

3.47 

256  x  256 

1.168  x  10~8 

2.92 

4.424  x  10~8 

2.83 

Table  6:  Accuracy  results  for  the  algorithm  described  in  section  4. 


presented  in  [9]  and  the  algorithm  described  in  section  4.  Figure  17  depicts 
the  interface  evolution  at  several  times  (top),  and  the  cross  section  of  the 
temperature  at  initial  (bottom  left)  and  final  (bottom  right)  times. 


Number  of  Points 

Ll  —  error 

order 

L°°  —  error 

order 

16  x  16 

8.047  x  10~4 

— 

2.709  x  10-3 

— 

32  x  32 

4.384  x  10-4 

0.876 

1.528  x  10~3 

0.826 

64  x  64 

1.874  x  10~4 

1.23 

9.724  x  10~4 

0.652 

128  x  128 

9.606  x  10“5 

0.964 

5.500  x  10"4 

0.822 

256  x  256 

4.723  x  10”5 

1.02 

2.822  x  10“4 

0.963 

Table  7:  Accuracy  results  for  the  algorithm  described  in  [9]. 


5  Conclusions  and  Future  Work 

We  have  proposed  a  simple  finite  difference  algorithm  for  obtaining  fourth 
order  accurate  solutions  for  the  Laplace  equation  on  arbitrary  domains.  We 
also  designed  a  fourth  order  accurate  scheme  for  the  heat  equation  with 
Dirichlet  boundary  conditions  on  an  irregular  domain.  In  the  case  of  the 
heat  equation,  we  utilize  an  implicit  time  discretization  to  overcome  the 
stringent  time  step  restrictions  associated  with  explicit  schemes.  We  then 
constructed  a  third  order  accurate  method  for  the  Stefan  problem.  We 
presented  multidimensional  results  to  demonstrate  the  accuracy  in  the  L°°- 
norrn.  Notably,  we  remark  that  in  two  spatial  dimensions,  one  can  obtain  6 
digits  of  accuracy  (for  the  Laplace  and  heat  equations,  5  digits  of  accuracy 
for  a  Stefan  problem)  on  very  coarse  grids  of  32  grid  nodes  in  each  spatial 
dimension.  Therefore,  even  though  the  discretization  yields  a  nonsymmetric 
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linear  system,  the  ability  of  this  algorithm  to  perform  well  on  a  very  coarse 
grid  makes  it  exceptionally  efficient.  Future  work  on  this  subject  will  include 
the  use  of  adaptive  mesh  refinement  techniques. 

The  Gibbs-Thonrpson  interface  condition  can  be  used  to  account  for  the 
deviation  of  the  interface  temperature  from  equilibrium  with  Tj  =  —ecn  — 
evVn,  where  k  is  the  curvature  of  the  front,  ec  the  surface  tension  coefficient, 
ev  the  molecular  kinetic  coefficient  and  Vn  the  interface  velocity.  Anisotropic 
surface  tension  effects  can  be  included  in  the  coefficient  ec.  For  example,  in 
two  spatial  dimensions,  one  can  take  ec  =  (1  —  15ecos(4a))  where  e  is  the 
anisotropy  strength,  and  a  is  the  angle  between  the  normal  at  the  interface 
and  the  x-axis.  In  [9],  the  Gibbs-Thonrpson  relation  was  computed  at  every 
grid  point  neighboring  the  interface  and  then  linearly  interpolated  to  the 
front.  In  that  case,  the  computation  of  the  interface  curvature  was  performed 
with  standard  second  order  accurate  central  differencing.  Such  computations 
for  the  curvature,  which  do  not  hinder  the  first  order  accuracy  of  the  method 
presented  in  [9],  cannot  be  used  in  this  present  work  without  lowering  the 
third  order  accuracy  of  our  method.  As  a  consequence,  more  research  on 
high  order  accurate  curvature  discretizations  must  first  be  performed  before 
considering  the  more  general  Gibbs-Thonrpson  case.  See,  for  example,  the 
work  of  [32]. 
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Figure  1:  Definition  of  the  ghost  cells  with  linear  extrapolation.  First,  we 
construct  a  linear  interpolant  T(x)  =  ax  +  b  of  T  such  that  T( 0)  =  T*  and 
T(9Ax)  =  Tj.  Then  we  define  T^1  =  T(Ax)  (likewise,  T^_2  =  T(2Ax)). 
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Figure  2:  Error  analysis  in  the  L°°-norm  for  the  one  dimensional  Laplace 
equation.  The  open  symbols  represent  the  errors  versus  the  number  of  grid 
nodes  on  a  loglog  scale,  and  the  solid  lines  are  the  least  square  fits  with 
slopes  -1.03,  -2.07,  -2.91  and  -4.10  in  the  case  where  the  ghost  cells  are 
defined  by  constant,  linear,  quadratic  and  cubic  extrapolation,  respectively. 


22 


3 

2.5 
2 

1.5 
1 

0.5 
0 

-0.5 
-1 
-1.5  — 

/ 

1 


Figure  3:  Solution  of  the  Laplace  equation  on  an  irregular  domain  in  two 
spatial  dimensions.  The  exact  solution  is  T  =  sin(7rx)  +  sin(7ry)  +  cos(7 rx)  + 
cos(7 ry)  +  x6  +  y6  in  .  The  grid  size  is  64  x  64,  and  the  ghost  cells  are 
defined  using  cubic  extrapolation. 
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Figure  4:  Error  analysis  in  the  L°°-norm  for  the  two  dimensional  Laplace 
equation.  The  open  symbols  represent  the  errors  versus  the  number  of  grid 
nodes  on  a  loglog  scale,  and  the  solid  lines  are  the  least  square  fits  with  slopes 
-.85,  -1.94,  -2.94  and  -3.96  in  the  case  where  the  ghost  cells  are  defined  by 
constant,  linear,  quadratic  and  cubic  extrapolation,  respectively.  Note  that 
on  a  32  x  32  mesh,  the  error  for  the  quadratic  extrapolation  is  the  same 
as  that  for  the  cubic  extrapolation.  This  is  an  example  where  there  was 
not  enough  points  within  the  domain  to  construct  a  cubic  interpolant  and 
the  algorithm  is  temporarily  forced  to  use  a  quadratic  interpolant.  More¬ 
over,  this  is  exacerbated  by  our  use  of  the  L°°-norm.  The  L 1  error  is  more 
forgiving,  since  it  is  based  on  averaging. 
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Figure  5:  Error  analysis  in  the  L°°-norm  for  the  one  dimensional  heat  equa¬ 
tion.  The  open  symbols  represent  the  error  versus  the  number  of  grid  nodes 
on  a  loglog  scale,  and  the  solid  line  depicts  the  least  square  fit  with  slope 
-4.14  in  the  case  where  the  ghost  cells  are  defined  by  cubic  extrapolation. 
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Figure  6:  Solution  of  the  heat  equation  on  an  irregular  domain  in  two  spatial 
dimensions.  The  exact  solution  is  T  =  e~2t  sin(x)  sin(y)  in  .  The  grid 
size  is  64  x  64,  and  the  ghost  cells  are  defined  by  cubic  extrapolation. 
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Figure  7:  Error  analysis  in  the  L°°-norm  for  the  two  dimensional  heat  equa¬ 
tion.  The  open  symbols  represent  the  error  versus  the  number  of  grid  nodes 
on  a  loglog  scale,  and  the  solid  line  depicts  the  least  square  fit  with  slope 
-3.94  in  the  case  where  the  ghost  cells  are  defined  by  cubic  extrapolation. 
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Figure  8:  Interface  at  time  tn  (dashed  line)  and  tn+1  (solid  line).  The  dark 
point  with  a  question  mark  represents  a  grid  node  that  is  swept  over  by  the 
interface  between  two  consecutive  time  steps,  and  where  a  valid  value  of  Tn 
needs  to  be  extrapolated  in  a  non-Cartesian  normal  direction  in  order  to 
evaluate  fn  =  Tn  +  ^An(A)Tn  in  equation  (7). 
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Figure  9:  Top:  Plots  of  several  isocontours  of  the  solution  extrapolated  in  the 
normal  direction  from  inside  the  disk  to  the  outside  (cubic  extrapolation). 
Bottom:  Exact  solution  for  the  sake  of  comparison. 
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Figure  10:  Error  analysis  in  the  L°°-norm  for  the  two  dimensional  definition 
of  the  ghost  cells  defined  by  cubic  extrapolation  in  the  normal  direction  as 
proposed  in  [1],  The  open  symbols  are  the  errors  in  the  L°°-norm,  and  the 
solid  line  is  a  least  square  fit  with  slope  —4.38. 
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Figure  11:  Error  analysis  on  the  effect  of  perturbing  the  velocity  by  an 
O(Ax)  (circles),  0( Ax2)  (squares),  and  0(Ax3)  (triangles)  amount  in  the 
case  of  the  one  dimensional  Stefan  problem.  The  symbols  represent  the 
numerical  solution  on  a  loglog  scale,  and  the  solid  lines  are  the  least  square 
fits  with  slopes  -.99,  -2.03  and  -3.05. 
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Figure  12:  Error  analysis  in  the  L 1  (circles)  and  L°°  (triangles)  norms  for 
the  two  dimensional  reinitialization  equation  (10).  We  use  a  fifth  order 
accurate  WENO  discretization  in  space  and  a  third  order  accurate  TVD 
Runge-Kutta  discretization  in  time.  The  initial  data  is  4>  =  ep~2'313  —  1, 
with  p  =  \J x2  +  y2.  The  symbols  represent  the  errors  on  a  loglog  scale,  and 
the  solid  lines  are  the  least  square  fits  with  slopes  -2.04,  -1.88,  respectively. 
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Figure  13:  Error  analysis  in  the  L°°-norm  for  the  one  dimensional  Stefan 
problem  on  a  loglog  scale.  The  ghost  cells  are  defined  via  cubic  extrapola¬ 
tion,  and  we  use  the  Crank-Nicholson  scheme  with  At  ~  Ax2.  The  stars 
depict  the  errors  when  the  exact  velocity  is  given,  and  the  triangles  illus¬ 
trate  the  errors  when  the  velocity  is  computed.  The  solid  lines  are  the  least 
square  fits  with  slopes  -4.03  and  -3.10. 
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Figure  14:  Error  analysis  in  the  L°°-norm  for  the  one  dimensional  Frank 
spheres  solution  (the  interface  velocity  is  not  constant  in  time).  The  ghost 
cells  are  defined  by  cubic  extrapolation,  and  we  use  the  Crank-Nicholson 
scheme  with  At  ~  As3/2.  The  symbols  represent  the  numerical  solution  on 
a  loglog  scale,  and  the  solid  line  is  the  least  square  fit  with  slope  -2.18. 
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Figure  15:  Error  analysis  in  the  L°°-norm  for  the  one  dimensional  Frank 
spheres  solution  (the  interface  velocity  is  not  constant  in  time).  The  ghost 
cells  are  defined  via  cubic  extrapolation,  and  we  use  the  Crank-Nicholson 
scheme  with  At  ~  Ax3/2.  The  time  discretization  involves  iterating  on  the 
velocity  as  described  in  section  4.3.  The  symbols  represent  the  numerical 
solution  on  a  loglog  scale,  and  the  solid  line  is  the  least  square  fit  with  slope 
-3.02. 


35 


Figure  16:  Interface  evolution  in  the  case  of  the  two  dimensional  Frank 
spheres  solution. 


Figure  17:  Cross  section  of  the  two  dimensional  Frank  spheres  solution  at 
initial  =  1  (left)  and  ffinal  =  2.89  (right). 
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Figure  18:  Error  analysis  in  the  L°°-norm  for  the  two  dimensional  Frank 
spheres  solution.  The  triangles  illustrate  the  accuracy  obtained  with  the 
scheme  presented  in  [9],  and  the  circles  represent  the  accuracy  obtained 
with  the  algorithm  presented  in  section  4.  The  open  symbols  are  the  errors 
in  the  maximum  norm,  and  the  solid  lines  are  least  square  fits  with  slopes 
-.80  and  -3.07,  respectively. 
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